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SUMMARY 


A three-dimensional, finite element computer code was developed to analyze ejector and axisymmet- 
ric fluted mixer systems whose flow fields are not significantly influenced by streamwise diffusion ef- 
fects. A two equation turbulence model was used to make comparisons between theory and data for 
variqus flow fields which are components of the ejector system, i.e., (a) turbulent boundary layer in a 
duct, (b) rectangular nozzle (free jet), (c) axisymmetric nozzle (free jet), (d) hypermixing nozzle 
(free jet), and (e) plane wall jet. Likewise, comparisons of the code with analytical results and/or 
other numerical solutions were made for components of the axisymmetric fluted mixer system. These 
included: (a) developing pipe flow, (b) developing flow in an annular pipe, (c) developing flow in an 
axisymmetric pipe with conical center body and no fluting, and (d) developing fluted pipe flow. 

Finally, two demonstration cases are presented which show the code’s abihty to analyze both the ejec- 
tor and axisymmetric fluted mixers. 


INTRODUCTION 

Solution of the compressible three-dimensional boundary layer equations for flows over arbitrarily 
shaped external and internal configurations is a current design and analysis requirement for numerous 
important systems applications. It is only in the last few years that efficient numerical techniques 
coupled with large high speed computers have become available thereby making solution of these 
problems feasible, e.g., see the Proceedings of the Conference on Aerodynamic Analysis Requiring Ad- 
vanced Computers, 1975 (ref. 1). Three-dimensional boundary layer solutions have been used to aid in 
understanding and/or designing both internal and external flow environments. For example, three- 
dimensional boundary layer solutions have been used to characterize mixing in scramjet combustors and 
chemical laser optical cavities, Zelazny et al. (refs. 2 and 3). 

The recent development of hypermixing nozzles (refs 4 and 5) for ejectors has made the application 
of ejectors to VTOL aircraft appear attractive. For example, ejectors are being used in the wing and 
canard of supersonic aircraft to provide the thrust augmentation for VTOL capability (ref. 6). Also, 
fluted nozzle mixers are being used in various systems (turbine exhaust, etc.) to conserve fuel by mixing 
a hot core flow with a cold outer flow, thus getting more thrust from the resulting homogeneous exit 
mixture. 

The overall objective of the present investigation, therefore, was to develop a computer program 
having the demonstrated capability of solving three-dimensional, compressible, viscous, mixing flow 
problems such as those characterizing the ejector system shown in figure 1 , and the axisymmetric 
fluted mixer shown in figure 2. This goal was accomplished by extending and modifying the existing 
COMOC (Computational Mechanics of Continua) finite element program developed in earlier studies 
(refs. 2, 7, 8, and 9) to the ejector and axisymmetric fluted mixer system classes. 

The report is subdivided into six sections plus the Appendices and References. The nomenclature 
used immediately follows this introduction. 

The governing conservation equations and models for the effective diffusion coefficients (laminar 
plus turbulent) are presented in the METHOD OF ANALYSIS. Also presented in this section are a 
description of the generalized boundary conditions and the coordinate transformation used to analyze 
geometries whose cross section varies with distance downstream. The finite element solution algorithm 
used to solve the governing equations is described in detail in reference 2. 
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The ejector system is comprised of several component parts (turbulent wall jet, hypermixing nozzle, 
etc.). Experiinental data characterizing these components were used to compare against theoretical 
predictions obtained using the COMOC computer code. A description of these experimental results 
and the details of the theoretical predictions are presented in the NUMERICAL RESULTS section. 
Also presented in this section are comparisons of the COMOC code with other numerical codes that 
model certain aspects of the axisymmetric fluted mixer system. Key results and recommendations for 
future efforts are presented in CONCLUDING REMARKS. Finally, APPENDICES A and B are in- 
cluded wherein the required input and output control for the COMOC code is described in a card-by- 
card sequence. 
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NOMENCLATURE 

Emperical constant = 23.5 

Emperical constants = 1.44, 1.92 

Skin friction coefficient 

Specific heat of ith species 

Dissipation rate of turbulent kinetic energy 

Static enthalpy of ith species 

Heat of formation of ith species 

Total enthalpy 

Turbulent kinetic energy (TKE) 

Mixing length 
Length scale 

Emperical constants = 1.0, 1.3 
Effective Prandtl number 

Reynold’s number 
Effective Schmidt number 

Radial coordinate 
Static temperature 
Total temperature 
Axial velocity in x-direction 
Transverse velocity in y-direction 
Transverse velocity in r-direction 
Transverse velocity in 6 direction 
Transverse velocity in z-direction 
Molecular weight of ith species 
Rectangular coordinates 
Mass fraction of ith species 
Boundary layer tliickness 
Displacement thickness 
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Edge of constant shear stress region 
Kinematic eddy viscosity 
Azimuthal coordinate, momentum thickness 
Von Karman constant = 0.435 
Emperical constant = 0.09 
Laminar viscosity 
Effective viscosity - ix + p e 
Kinematic viscosity 
Density 
Shear stress 

Von Driest Damping Coefficient 
Coordinate system in transformed domain 
Evaluated at edge of boundary layer 
Turbulent value 
Reference value 

Nondimensional quantity in wall shear stress formulation 
Effective quantity, equal to sum of laminar and turbulent value 
Nondimensional quantity in annular pipe formulation 
Quantity evaluated at y = 6 
Differentiation with respect to ^ 

METHOD OF ANALYSIS 

Many researchers are now giving attention to numerical solution of three-dimensional parabolic and/ 
or boundary-region flow fields. Most procedures employ a finite difference solution algorithm for 
variously combined forms of the continuity, momentum and energy equations. Several researchers 
have obtained solutions for the three-dimensional boundary-region flow of single-species fluids. A brief 
description of several methods (refs. 10 to 17) that have been employed is given in reference 2. 

The system of partial differential equations governing such three-dimensional, confined unidirec- 
tional flows of a compressible, reacting fluid is obtained as an approximation to the full three-dimen- 
sional Navier-Stokes equations. This approximation, referred to as the “parabolic Navier-Stokes 
equations,” describes steady, confined three-dimensional flows wherein (1) a predominant flow direc- 
tion is discernible; (2) diffusion of mass, momentum, and energy in the predominant flow direction is 
negligible compared with convection; and (3) no disturbances are propagated upstream, e.g., recircula- 
tion is not considered. 

Solution of the three-dimensional parabolic Navier-Stokes equations is obtained wherein the stream- 
wise pressure variation is computed using the approach presented in the subsection Pressure Variation 
in Ducted Flows. This version of the computer code which computes an axial pressure gradient for 
confined flows is referred to as COMOC-3DPNS. If a pressure distribution is known a priori then that 
section of the code which computes streamwise pressure variations may be bypassed. The resulting 
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equations are herein referred to as the three-dimensional boundary region equations and represent a 
subset of the 3DPNS system. This variant of the code is referred to as COMOC-3DBR. Details of the 
numerical solution algorithm used to solve the differential equation system may be found in references 
2 and 8. 


Governing Equations 

The parabolic Navier-Stokes equation system for a multiple-species, compressible, nonreacting flow 
takes the form: ' 
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Species Continuity 


‘J^i ‘’’^i a / ‘'MX a / or 

nil + pv + PW = — / 1 + — / 

ax ay dz ■ ay Ij^e dy f . dz I j^je jsj dz 

\ Sc Re f \ Sc 


For parabolic flow in the x-direction, p in equation (2) is the mean value over the local (y, z) cross 
section and the diffusion terms in the x-direction have been neglected implying that recirculation zones 
are negligible in the flow fields of interest. The variables appearing in equations (1) to (6) are non- 
dimensionalized with respect to Poo, and a length scale L. The Reynolds number 

effective Prandtl number N® , and effective Schmidt number Nc„ are defined for a combination of 
laminar and turbulent contribuiions as, 

NRe = °° - (7a) 


P + 


jyje Npj. (Npj.) 

Pr T 


a: _££_ (7e 

N"sc Nsc 

In equation (7), p is the laminar viscosity, e is the kinematic eddy viscosity, and the subscript T de- 
notes a turbulent parameter. The stagnation enthalpy is defined in terms of species static enthalpies 


H = 2 hjYj -H — (u^ + v^ + wM (8) 

i ^ 

The static enthalpy includes the heat of formation h^ j of the ith species in its definition as 

hi = / Cp,idt + (9) 

^o 

An equation of state is obtained by assuming perfect-gas behavior for each species and using Dalton’s 
Law, 




where R is the universal gas constant and W: is the molecular weight of the ith species. 


5 


Effective Diffusion Coefficient Models 


Closure of the governing equations requires introducing relationships to define the effective viscosity, 
and the turbulent Prandtl, and Schmidt numbers. In this investigation it was assumed that mass and 
thermal energy, both scalar quantities, diffuse at the same rate. Hence, the turbulent Lewis number, 
(NLe).p! is equal to unity. The validity of this assumption is supported by numerous experimental 

studies and is generally used in most turbulence modeling efforts. This assumption requires (Npj.)^ = 


(Nsc).j,- It is then only necessary to specify (Np,.)^. Generally, it is sufficient to define (Npj.)^ equal 
to a constant in the range from 0.7 to 1 .0. For this study, the value of (Npj.)^ = 0.7 is used. 


Four different approaches to modeling e are considered: (1) e equal to a constant (usually some 
multiple of the laminar viscosity, p), (2) mixing length theory (MLT), (3) a one-equation turbulence 
model coupled with an empirical turbulent length scale relation, and (4) the two-equation turbulence 
model reported by Launder and Spalding (ref. 18). 


Mixing Length Theory - The two-dimensional MLT in the (x, y) coordinates, e is given by: 
3u 


- 


3y 


( 11 ) 


The mixing length is defined as 


Kyco (0 S y s X5 / «) 

C = 

XSo) (y > X 5 /k) 

where 

K = 0.435 

X = 0.09 

8 = boundary-layer thickness 

The Van Driest damping coefficient is 
y * 

CJ = 1 - exp (—-) 

A 

where: 


u* = friction velocity, >/ 

T = skin friction 

p = density at wall 


( 12 ) 


(13) 
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V 


kinematic viscosity 


A = 23.5 

One (k) and Two (k&d) Equations Turbulence Models - Dimensional analysis shows that the kinematic 
viscosity may be expressed as the product of a characteristic velocity (V^) and the length scale 
henc6 define 


e = (14a) 

where is an empirical constant and the relationship between and £ from MLT is given by 
Launder and Spalding (ref. 1 8) as 

£d = (14b) 

The two models both assume that the turbulence velocity is given by 

Vt = k*/^ (15) 

where k is the turbulent kinetic energy 


k = 



+ 


v'^ + w'2 
2 


(16) 


The one equation model assumes an empirical relation for £^j, while the two equation model 
expresses in terms of k and the dissipation rate of turbulence, d as 


= k3/2/d (17) 

The turbulent kinetic energy and dissipation are computed from 
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(18) 


(19) 


where it is assumed that the production of turbulent kinetic energy is due to the u' v' component of 
shear stress, i.e., 0u/0z is negligible compared to 0u/0y. This assumption is reasonable for the flow 
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regimes examined herein using the k and d model. The emperical constants used in equations (14), 
(18), and (19) as suggested in reference 18 are as follows: C, = 0.09, cf = 1.44, = 1.92, Nu = 1.0, 

and N^i - 1.3. 

Pressure Variation in Ducted Flow 


For internal flows, characterized by boundary-layer thicknesses which are small in comparison with 
the overall internal duct dimension, the axial pressure distribution can be accurately approximated by 
inviscid flow solutions. However, when the flow is confined in a duct whose lateral dimension is not 
large with respect to the boundary layer thickness, this approach is invalid. Here, boundary-layer 
development directly influences the pressure distribution within the duct, and an axial pressure 
gradient is induced by viscous effects. For these flows, a quasi-one-dimensional integral treatments 
of equations (3) and (5) has been suggested wherein for steady flows, equations (3) and (5) are inte- 
grated across the duct transverse dimensions to obtain a relation for the mean value of axial pressure 
gradient (three-dimensional parabolic Navier-Stokes variant). Details of the method are described in 
references 17 and 19. 

Generalized Boundary Conditions 

The finite element solution algorithm considers the following generalized boundary condition for 
the variable q. 

dq 

ai9 + a2 = 0 (20) 

where a j , 32, 33 are user specified constants and dq/dn is the local outward pointing normal deriva- 
tive. 


For turbulent flows near walls where the velocity gradient becomes very large, it was found useful 
to replace the no-slip boundary condition for the axial velocity. In its place, the following boundary 
condition on the velocity gradient at the wall was used. 


du \ _ 
''''wall ' " 


( 21 ) 


where is the local value of the wall shear stress and is computed from the following relation given 
by Patankar and Spalding (ref. 19). 
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.. ( 22 ) 


where k is von Karman’s constant (used in defining the mixing length) R^ = 


k 2 R, R = — The 

t^l 


subscript 1 corresponds to evaluation at a distance 6 from the wall (6 is defined to be the edge of 
the constant shear stress region). The effect of pressure gradient on the shear stress is included in 
equation (22) by the parameter Fp given by 
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1 - 4F R 
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dp 
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(23) 


Variable Geometry Coordinate Transformation - Rectangular 


By an appropriate transformation of variables (essentially coordinate stretching) the cross section 
in the (y, z) plane may be allowed to grow as a specified function of x. This proves useful in com- 
puting internal flows (wall jets, etc.) when the rate of growth of the solution field can be approximated 
a priori. It is also essential in computing confined flows in diffusers or other variable area ducts. It 
should be noted that the finite element grid is established at the initial cross-sectional geometry (what- 
ever its shape) and then the finite elements in the cross-section expand and/or contract as the flow is 


njiiLi'wctj uwWiistreaiii. At the present timC, the COMOC code can handle constant area flows with 


any cross-sectional geometry. However, the variable cross section option presently is limited to 
rectangular or axisymmetric geometries. 


The coordinate transform for the ejector system is defined as follows. Consider a flow in a domain 
bounded in the y direction by walls with a specified geometry y = f) (x) and y = f 2 (x) and bounded 
in the z direction by walls z = g| (z) and z = g 2 (x) (see Figure 3). The transformations from the 
(x, y, z) space to the (|, T),f) space are given by; 


? = X 
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^ f 2 (x) - fj (x) 


z - gi (x) 


r - 
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Then, for example, equation (2) in the transformed plane becomes 
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(I) 


§2 (^)-gl (0 


F3 a) = 


G 3 (i) = 


f2(?)-fi(l) g2(?)-gl(l) 

Variable Geometry Coordinate Transformation - Axisymmetric 

. The coordinate transform for the axisymmetric fluted mixer system is performed using the same 
approach used for rectangular geometries. Consider a flow in an axisymmetric domain bounded in 
the radial direction by walls with a specified geometry (see figure 4). The transformation from the 
(x, y, z) to the (?, tj, f) domain is then found as follows. All points in the (y, z) plane can be specified 
in polar coordinates as (see figure 4). 


y = rsin0 
z = rcos0 

Likewise, the corresponding point in the corresponding (j?, f) plane is given by; 
77 = rsind 
f = rcos0 
where 
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where; 
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NUMERICAL RESULTS 

A good test of the accuracy of a numerical model is a direct comparison of theoretical predictions 
with experimental data. To this end, predictions for nine different flow conditions were obtained, 
table II. These particular test cases were chosen for one or both of the following reasons: the data 
were well documented in that turbulence initial conditions for the TKE equations were either speci- 
fied or obtainable using algebraic relations from other given data; a given case represents one of the 
possible component parts of the more complex ejector or axisymmetric fluted mixer system. Where 
applicable, more than one of the turbulence models described above was used. 

Test Cases 

Case 1. Turbulent Boundary Layer in a Duct - The data for this case were selected from the Proceed- 
ings of the Stanford conference on the Computation of Turbulent Boundary Layers (ref. 2 1 ) and pro- 
vide a well defined description of incompressible boundary layer development. The specific case 
(Number 2400 of ref. 21) represents an initial equilibrium boundary layer in a moderate positive pres- 
sure gradient (specified) where the pressure gradient is abruptly decreased to zero and the flow relaxes 
towards a new equilibrium condition. Since the pressure distribution is given, the three-dimensional 
boundary region variant of the COMOC computer code was used. The initial velocity profile was used 
to generate initial conditions for e'and using the MLT model. The initial k&d profiles were then 
obtained algebraically from equations (14), (1 5) and (17). The gradient boundary condition relating 
the velocity gradient at the wall to the local value of the wall shear stress was used. The transverse 
velocity component v is obtained by numerically integrating the two-dimensional form of the con- 
tinuity equation. 
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' Data available for comparison include the axial (x) development of the velocity profile, skin friction 
coefficient C.f/2, displacement thickness 5*, and momentum thickness d, where: 


Cf/2 = 


w 

pu| 


5* - (1 ._^)dy 


5 pu u 

e = /o — (1 - -)dy 

PgUe u. 


(33) 


(34) 

(35) 


By definition, 6 is the local boundary layer thickness defined as the distance y from the wall to the 
point in the boundary layer profile where u/ug equal 0.99 and Ug,Pe are the inviscid values at the 
outer edge of the boundary layer. 

Results were obtained using the MLT model and the two-equation TKE model. Figure 5 shows the 
skin friction Cf/2, 5* and 0, while figure 6 depicts the final velocity profile at 241 .3 cm downstream. 
As can be seen, both models accurately predict (within 5%) both the final shape of the velocity pro- 
file as well as the parameters describing its axial development. Figure 7 shows a comparison of two 
runs made, both of which used the MLT model. However, one run was made using the gradient bound- 
ary condition. The finite element discretization for both cases was the same. As can be seen, the pre- 
diction of 5* and 6 are within 5% of the data for both cases. However, the skin friction coefficient 
Cf/2 drops off very rapidly using the no-slip boundary condition (at a distance of 132.1 cm, the pre- 
dicted value is less than half the measured value). Better agreement could be obtained using the no 
slip wall only by adding several more nodes in the near wall region. This would greatly increase the 
amount of computer time required since the integration step size stability is controlled by the mini- 
mum spacing between the nodes in the transverse direction. This emphasizes the benefits obtained in 
using the gradient boundary condition. 

Case 2. Rectangular Nozzle (Two-Dimensional Flow) - The data for this case is given by Bradbury 
(ref. 22). He obtained self-similar velocity and TKE profiles by empirically curve fitting his experi- 
mental data. He gives the following results: 


HJLHi = f(^) = exp [-0.6749tj2 ( 1 + 0.0269tj'‘ ) ] (36) 

= g(7?) (37) 

“o 

where 

T? = y/5 

Uq ~ centerline velocity inside jet 

u, ~ freestream velocity 
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5 ~ y location inside jet where f(Tj) = 1/2 

g(r]) ~ given in graphic form 

The computations were started at x/D equal to 20 (where D is the nozzle exit width in the y-direction). 
Initial profiles for u and k were obtained from the empirical curves. The initial value of 2(j was ob- 
tained^from the empirical relation: 

2d = Iy 2 -yil (38) 

where 

yi = point in jet profile where f(r?) =0.1 

Y 2 = point in jet profile where f(7j) = 0.9 

Then the initial d profile may be obtained using equation (17). Again, the transverse velocity is com- 
puted by numerically integrating the two-dimensional form of the continuity equation. 

Results were obtained using the two-equation TKE model and also using the TKE equation ( 1 8) 
combined with equation (38), from which 2d is computed at each x station based on the local shape 
of the velocity profile. Figure 8 shows results using both models for the centerline velocity decay and 
the half-jet thickness growth compared with the experimental data. The two-equation model gives 
excellent agreement for Uq/ui and overpredicts S/D by less than 5% at x/D equal to 70. The one-equa- 
tion model overpredicts Ug/ui by 6% and underpredicts 6/D by 7% at x/D equal 70. However, the 
true test of both models is in how well they compare with the turbulence data. 

Figure 9 gives the results of the centerline turbulence kinetic energy (k^) scaled by Uq for both 
models. The two-equation model shows an initial drop in this value then a steady increase so that 
at x/D equal 70 the value is 8% higher than the experimental value. The one-equation model shows 
an initial drop and then the value remains constant at 10% lower than the experimental value. The 
discrepancy in both models may be due to the values of 0.1 and 0.9 used to compute 2(j in equation 
(38). These values are poorly defined and further refinement could give better agreement for both 
models since these values are used to generate an initial d profile for the two-equation model. The 
shape of the k profile for both models at x/D equal 70 is shown in figure 10. The profiles have been 
scaled with respect to k^. so that the magnitude differences depicted in figure 9 are eliminated. Both 
models accurately predict the shape of the profile (to within 1 0%) as well as the location of the maxi- 
mum turbulence kinetic energy. 

Case 3. Axisymmetric Nozzle - This case provides a good check on the three-dimensional solution 
capability of the analysis. Although the flow field may be characterized using only two spatial co- 
ordinates (the radius r and axis x), solutions were obtained using a rectangular coordinate system 
(x,y,z). The specific case considered is an incompressible jet into a stationary medium (quiescent 
jet). Calculations were started at 106.7 cm downstream of an 8.89 cm diameter jet (x/D = 12) where 
the flow may be described using the similarity solutions reported by Schlichting (ref. 23). The sim- 
ilarity solution is vahd for a constant effective viscosity e. Therefore, solutions were obtained by 
integrating the three velocity components using a constant value of e and also the two-equation TKE 
model. Initial values for the two-equation model were obtained algebraically by using the known 
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value of € and by estimating the initial value of from the velocity profile. Figure 1 1 shows the 
predictions for both models compared to the similarity solution. 

Case 4. Hypermixing Nozzle - Figure 12 shows the hypermixing nozzle being considered. Data are 
available from Eastlake (ref. 24) giving velocity profiles in the y, z plane at distances 3.0 and 5.0 in. 
downstream from the nozzle With the exception of the jet exit speed and nozzle geometry, no initial 
condition data were available for this case. The procedure used to establish initial conditions was to 
first assume a constant value of e for the turbulence model. This value was varied until reasonable 
agreement was obtained with the experimental data for the maximum velocity decay. It was assumed 
(arbitrarily) that the turbulence kinetic energy in the hypermixing nozzle exit plane was 1 .5% and 
equations (14), (15), and (17) were used to compute the initial turbulence dissipation. Having defined 
the starting conditions, calculations were then made to predict the flow field development down- 
stream using the two-equation turbulence model. Both magnitude and shape of the velocity profiles 
in transverse direction were reasonably predicted, figure 13, with a maximum difference between data 
and theory of 5% for the peak velocity. Of course, these calculations only demonstrate the capability 
of the code to analyze the hypermixing configuration rather than validate the turbulence models or 
manner in which the initial conditions were obtained. Analysis of turbulence data for this configura- 
tion would be required to address the question relating to the accuracy of the turbulence models and 
initial conditions. Such an effort was beyond the scope of this study. 

Case 5. Plane Wall Jet - This flow provides an excellent check of the turbulence model being used. 
Consider a two-dimensional flow where the turbulent component of the shear stress is assumed to be 
given by: 


3u 

~ (39) 

An obvious implication of equation (39) is that the turbulent stress goes to zero at the maximum 
velocity location in the u(y) profile for a turbulent wall jet. It has been pointed out experimentally 
by Harris (ref. 25) that the turbulent stress does not go to zero at this point. Indeed, the turbu- 
lent stress is quite large and is caused by the asymmetry of the velocity profile about the maximum 
velocity point. 

The plane turbulent wall jet is shown in figure 14, along with a definition of the flow parameters. 
Initial profiles for u, k, and d were obtained as follows. Experimental data are available for self- 
similar velocity profiles (ref. 26), mean turbulent kinetic energy (ref. 27), and turbulent shear stress 
profiles (ref. 28). From these profiles and equations (14) and (17), initial profiles for u, k, and d at 
x/L = 25 were obtained and then numerically integrated to x/L = 100. The solution domain was 
allowed to grow with x using the variable geometry coordinate transformation. In this way, at x/L 
equal 100 there were still the same number of node points in the turbulent boundary layer portion 
of the jet as at the initial station. Both the no slip wall and the gradient boundary condition were 
used to provide another check on the accuracy of the gradient boundary condition. 

2 

Shown in figure 1 5a are numerical results for the axial skin friction Cf given by iT^jpUy com- 
pared with experimental data of reference 26. As can be seen, the results obtained using the 
gradient boundary condition remain within the spread of the data. It should be pointed out that 
both numerical curves were obtained using only four node points in the initial u profile between 
the wall and the maximum velocity. This explains why the no slip boundary condition results for 
the skin friction are over an order-of-magnitude smaller than the data at x/L equal 100, i.e., the 
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discretization was inadequate. Shown in figure 1 5b are the velocity profiles at x/L = 1 00 for both 
cases compared with the data of Schwartz and Cosart (ref. 28). Again, the gradient boundary con- 
dition results show the flow field remaining self-similar, while the no-shp boundary condition results 
showed marked deviation from the self-similar profile. Figures 16a and 16b show the maximum 
velocity decay and half-jet growth compared with the data of reference 26. Figure 1 7a shows the 
maximum turbulent kinetic energy (kj^ax) compared to the experimental self-similar value. This 
value increases initially to a value 9% larger at x/L equal 30 then decreases to a value only 1 % larger 
at x/L equal 100. Figure 17b also shows the turbulent kinetic energy profile at x/L equal 100. 

Case 6. Developing Pipe Flow - The experimental data for this case are taken from Nikuradse (ref. 
23) who shows the development of the velocity profile from its initial slug value to the fully 
developed parabolic profile. For this case, a two-dimensional finite difference code in polar co- 
ordinates was also used to make comparisons between the predicted secondary velocity decay com- 
puted using finite element and difference techniques. 

For axisymmetric flows, the continuity equation is rewritten in polar coordinates, assuming the 
circumferential velocity component (v^) is negligible, to give 

0 0 

— — (pru) + — (prVf) = 0 (40) 


(41) 


(42) 


The radial velocity Vj. is related to v and w by 
V = Vj sind 

W = Vj. COS0 

where 

d = tan"‘ (-] 


The method of solution used by the COMOC code for axisymmetric flow is to integrate the u 
momentum equation (2) downstream, then compute Vj. from the continuity equation (40) in polar 
coordinates, and finally obtain v and w from equation (41). For cases 6, 7, and 8, the assumption 
of zero circumferential velocity, vg, is exact since the flows are two dimensional. However, case 9 
is three dimensional and will determine how much accuracy is lost when making the assumption that 
V0 = 0. The COMOC code was run for both a 45° and 90° wedge of the circular cross section pipe 
(using angular symmetry) with this method of computing v and w. The finite element discretization 
for the 90° wedge is shown in figure 18. The solution was started at x equal zero assuming a slug 

2pua 

profile for u and a flow Reynolds number (Re = ) of 500 where u is the initial slug velocity 

M 

and a is the pipe radius. The results for the growth of the centerline velocity uq^ , compared with 

the finite difference solution, are shown in figure 19a. The decay of the maximum radial velocity 
Vj- is shown in figure 1 9b. Both COMOC runs were made with the same number of elements. 
Therefore, the 45° wedge had finer resolution. As can be seen, both the 45° and 90° wedge give 
excellent agreement with the 2D finite difference code. 


Case 7. Developing Flow in an Annular Pipe — The analytical solution for a fully developed flow in 
an annular pipe is given by: 


■ (“)(S 


(1 - rJ^) 1 nr* 


(1 -r*)" 


Inr* 


(43) 
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where 


u* 

r* 

r* 

0 

X* 

Re 

p* 


= u/u 
= r/a 
To /a 
= x/a 
2pau 
M 
P 


( 44 ) 


and u is the initial slug velocity, a is the outer pipe radius, and ro is the inner pipe radius. 


dp 

The pressure gradient — may be found in terms of the mean velocity u, Reynolds number Re, 
dx 


and ro as follows. 



u 


= /' r*u*dr*//‘ 

r? r* 


r*dr* 


(45) 


using equations (43) and (44) gives 


i* = 


Re dp* 
16(1 -rg^) dx* 



Inr* J 


Therefore, 

^ ^ - 16(1 -r*^) 

dx* Re [(1 -rf ) + (l -r*')/lnr*] 


(46) 


(47) 


The two dimensional finite differences and COMOC codes were run with Re = 500, u = 0.5468, and 


ro = 0.25. The fully developed value of dp*/dx* = -0.0938 was obtained with both codes. 


Figure 20 shows the axial velocity profiles for both codes at x = 71a compared to the analytical 
solution. The COMOC solution domain is a 90° wedge with 5 1 node points (73 elements). The 
decay of the maximum radial velocity (w along 6 equal 0° and v along 9 equal 90°) is shown in 
figure 2 1 . 

Case 8. Developing Flow in an Axisymmetric Pipe with Conical Center Body and no Fluting — This 
case can be viewed as a combination of cases 6 and 7. At the initial x station, the pipe has an 
annular shape. This annulus then conically decreases in size until at some downstream station the 
pipe cross section changes from annular to circular. See figure 29b for a side view. For this problem, 
the variable geometry option for axisymmetric flows in COMOC was used to account for the 
conical center. The 2D finite difference code was also modified to compute variable geometry flows 
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in the radial direction. Shown in figure 22 is a comparison between the two codes for the u velocity 
profiles at several stations. The conical center has an initial radius of ro/a = 0.25 at x/a = 0, and 
its radius decreases linearly to zero at x/a = 10.0. Figure 23 shows the comparison of the secondary 
velocity decay between the two codes. 

Case 9. Developing Fluted Pipe Flow — The configuration for this flow is shown in figure 24. Two 
flows, of differing velocity, are separated by a fluting. At x = 0, the fluting ends and the two 
flows then merge in the downstream x direction. This flow is three dimensional in nature and was 
chosen to test the validity of the assumption that vg is negligible when computing Vj. from the 
continuity equation (40) in polar coordinates. A three dimensional finite difference code used in 
reference 29 was used as a comparison to test this assumption. While this code was written 
for more elaborate flows, it can be used to solve for the flow in a straight pipe in cyUndrical co- 
ordinates. Thus, both Vj. and vg may be computed for the configuration of figure 24. The results 
of these computations are shown in figure 25 where axial growth of the centerline velocity, non- 
dimensionalize by a reference velocity u, for both codes is presented. Figure 26 shows the u 
velocity profiles for both codes at several locations. As can be seen, both codes predict the fully 
developed parabolic profile at a distance of 50 pipe radii downstream. Figure 27 shows the key 
results for this case. Here the maximum radial velocity Vj- for both codes is presented as well as 
the maximum angular velocity \q for the finite difference code. As can be seen, the agreement for 
Vj. is excellent (less than 8% difference), while the finite difference code predicts a value of vg no 
larger than 6% of Vp Thus, the assumption that \g is negligible in equation (40) has been justified 
for this flow configuration. 

EJECTOR ANALYSIS 

Application of the code to aid in optimization of ejector design or compare theory with ejector 
data was beyond the scope of this investigation. However, a demonstrator case was considered to 
illustrate how COMOC might be employed to support ejector design studies. Consider a typical 
ejector configuration (ref. 30) composed of free jets, wall jets, a constant area mixing section, and 
a diffuser, figure 28. Each component of this system was analyzed in the test cases described above. 
Here our objective was to demonstrate the capability of the model to analyze the integrated com- 
ponents and their interaction. The exit plane (y, z) consists of a bank of nozzles (one of which is 
depicted in figure 28) spaced 7.62 cm in the z-direction. Each 25.4 cm high nozzle consists of a 
hypermixing nozzle (which is broken up into eight 2.54 cm segments) with a 2.54 cm rectangular 
nozzle at the top and bottom. Boundary layer control wall jets (0. 1 55 cm) are placed along the 
bottom and top walls. The mixing section is 25.4 cm long in the x-direction, followed by a diffuser 
with a total angle of 1 9.6 degrees. 

The desired information for this system is the thrust augmentation obtainable for various exit 
plane conditions (velocities, static pressure, temperature, etc.). The procedure used is as follows; 

(1 ) Exit plane conditions for the nozzle bank and wall jet are assumed to be given. 

(2) A guess is made for the value of the secondary (entrained) velocity. 

(3) The flow is then computed through the mixing section to the diffuser exit. 

(4) If the exit static pressure is correctly predicted to be ambient, then the thrust augmentation, 
skewness, etc. are computed for the system. If the exit pressure is not correctly predicted to be 
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ambient, then a new value for the secx)ndary velocity is selected and the calculation repeated. An 
alternate approach would be to calculate the length required to give ambient static pressure at the 
exit rather than fix the length of the diffuser section a priori. 

The initial conditions required at the nozzle exit plane are obtained in the following manner: 

( 1 ) The exit velocities (hypermixing, boundary layer control, and secondary) and the exit 
pressure were obtained using the one dimensional code developed by Salter for the same geometry 
(ref. 31). These values are: 

up ^ hypermixing velocity = 173.7 m/sec 

UfiL = boundary layer control velocity = 173.7 m/sec 

ug = secondary flow velocity = 51.8 m/sec 
Exit Static Pressure = 99626 N/m^ 

(2) It is assumed that the boundary layer control nozzle and the rectangular nozzle at the base 
of the hypermixing section possess fully developed turbulent wall jet and boundary layer profiles, 
respectively, at the jet exit. This assumption allows use of a relatively coarse grid spacing in the 
near-wall region at the nozzle exit plane and considerable savings in computer time due to the 
absence of steep velocity gradients. The initial velocity, turbulent kinetic energy (k), and dissipa- 
tion (d) profiles can then be determined from the approach previously used in cases 1 and 2 (see 
Numerical Results Section). 

(3) The initial k and d values for the hypermixing nozzle are obtained from the results of 
check case 4. The initial values of k and d for the present configuration were obtained by using 
the same nondimensional values of k/up and db/up from case 4, where b is the nozzle exit width. 

These values are: 

k/up = 1.46 x 10“^ 

db/up = 2.10x10*'* 

(4) The values for k and d in the secondary flow were nominally taken to be 10% of the hyper- 
mixing values. 

The ejector/diffuser system, figure 28, was analyzed by employing the COMOC restart option by 
making three separate runs each requiring 5 hours of IBM 360/65 computer time. Table I gives 
the initial and computed values of the velocities (u, v, w), turbulence quantities (k&d), etc. at x = 
0.0, 15.2, and 27.9 cm (see figure 28). 

AXISYMMETRIC FLUTED MIXER ANALYSIS 

Figure 29 illustrates the geometry for a 90° sector of an axisymmetric fluted mixer which was 
selected to establish the feasibility of using COMOC to characterize the three-dimensional flow 
field development. The outer fan flow, figure 29, is used to cool the inner core flow as it mixes 
downstream from the mixer exit plane. The radius of the center plug is assumed to decrease 
linearly from 12.7 cm at the exit plane to zero at 76.2 cm downstream. The radius of the outer 
cowl wall is assumed to be constant at 50.8 cm. The solution domain is continued an additional 
76.2 cm downstream from the end of the center plug. The initial conditions for this system are 
assumed to be: 
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Core Flow 


u = 260 m/sec 

To = 707°K 

Mach Number = 0.5 
P - = 41364 N/m^ 

Fan Flow 

u = 179 m/sec 

To = 333°K 

Mach Number “0.5 
P = 41364 N/m^ 

The inner and outer walls are assumed to be adiabatic. The inner wall axial velocity condition 
changes from no sUp on the centerbody wall to zero gradient (symmetry) on the centerline at the 
end of the center body. The turbulent diffusion rates were assumed to be equal to 1 00 times the 
laminar value rather than using a higher order turbulence model (represented a savings of 50% in 
computer time). This assumption is not required nor recommended when performing fluted mixer 
design studies. The results for the axial velocity and static temperature profiles along the radial • 
lines 0=0° and 0 = 22.5° at x = 0.0, 76.2 and 152.4 cm are shown in figures 30 and 31. 

CONCLUDING REMARKS 

A three dimensional computer code has been developed for analyzing ejector/diffuser and 
axisymmetric fluted mixer systems. Some specific results are; 

(1 ) A general three dimensional finite element model employing the two equation turbulence 
model was developed which characterizes flow field development not significantly influenced by 
streamwise diffusion effects. 

(2) This model was used to examine various flow fields and a comparison made between the 
model and data, i.e., (a) turbulent boundary layer in a duct, (b) rectangular nozzle (free jet), (c) 
axisymmetric nozzle (free jet), (d) hypermixing nozzle (free jet), (e) plane wall jet. 

(3) It was found that the initial conditions have a first order effect on the numerical results. 
Here a specific problem area is the initial turbulence parameter profiles required in the two 
equation turbulence model. 

(4) The CPU time for the nine test cases, the ejector system, and the axisymmetric fluted 
mixer presented above is presented in Table II. Also shown in this table are the number of in- 
tegrated variables, number of finite elements, and the number of integration steps for each case. 

It should be emphasized that the code can most efficiently be used to aid in providing design guid- 
ance to a specific device only after the theoretical predictions are compared to data representative 
of the device of interest. Two demonstration cases were analyzed to show the interaction of the 
different component parts for these devices. Detailed solutions do not come cheaply, i.e., 15 hours 
of computer time on the IBM 360/65 were required for the ejector demonstration case. 
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APPENDIX A 


COMOC DATA DECK PREPARATION 


A description of the data deck preparation for the 3DBR and 3DPNS versions of COMOC was re- 
ported in references 2 and 7, respectively. The 3DBR version is for flows where the pressure field is 
given and no area constraints are imposed upon the flow. The 3DPNS version reported in reference 7 
considered constant area ducted flows. In modifying the code to consider variable area, ducted injec- 
tor/diffuser systems and axisymmetric fluted nozzle mixer systems, several changes were required in 
the input and output subroutines, i.e., data deck preparation. For this reason, data decks for four 
different problems that typify these system classes have been chosen to demonstrate the required 
modifications. These decks and their pertinent features are: 


Four Test Cases 


Case 


Key Features 

I. 

Turbulent 

• 

Two-Equations (k&d) Turbulence Model 


Boundary 

• 

Gradient Boundary Condition 


Layer 

• 

Boundary Conditions for k&d at Node off the Wall 



• 

Boundary Layer Integral Parameter Output 

II. 

Diverging 

• 

Variable Geometry 


Channel 

Flow 

• 

Flow Separation 

III. 

Rectangular 

• 

Ejector Performance Parameters Output 


Nozzle 

Ejector 

• 

Constant Turbulent Viscosity Option 

IV. 

90° Sector of 

• 

Finite Element Setup and Variable Geometry for 


Axisymmetric, 


Nonrectangular Cross Section 


Conical 

• 

Generation of Cross Flow Velocities from Continuity 


Center, 
Fluted Mixer 
System 


Equations in Polar Coordinates 


Tables III-VI are a listing of the data decks used to run Cases I -IV. Tables VII-X are the output gen- 
erated by each of these decks. Details and a complete description of these tables are given below. 

DATA DECK DESCRIPTION 

Examining Tables III - VI, it is seen that the data deck may be conveniently divided into 27 separate 
sections. The purpose of each section, as well as the definition of terms as applied to each test case, 
are given below. 

Section 1. 

Starting in Column 1, FEBL, used to start execution of COMOC. 
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Section 2. 

Starting in Column 1, 3DPNS (Cases II, III, IV) for ducted flows where the pressure is computed 
internally or, 3DBR (Case I) in flows where the static pressure field is given. 

Section 3. 


Starting in Column 1, FENAME. This section consists of two parts for inputting integer and fixed 
point variables. For the first part, on the next card starting in Column 2, SNAMEOl , followed on 
succeeding cards with the input integer parameters. There is no designated column for these variables, 
as long as a comma follows immediately after each variable. A complete list of all the possible integer 
variables is given in Table XI, which was outputted for Case I. For the four test cases presented here, 
the required integer variables are as follows. The pertinent case (s) for each input variable are in paren- 


1 WXlIg l.Xl.V' V* W«. AllA MVy XA MAW < MA AM w > • . 


IGAS (I,II,III,IV) 
IFR (I,II,III,IV) 
KBNOTW (I) 
KBNOS (IV) 

KDUMP (I,II,III,IV) 
LG (IV) 


0 Default, isoenergetic flow with constant Cp 

1 General flow 

0 Default, equilibrium composition or complete reaction 

1 Frozen composition 

0 Default 

1 Use gradient boundary condition for variable 1 (Ul) 

0 Default 

n Number of initial nodes for Ul (KBNOS) and/or H (KBNOSH) 
whose boundary condition changes from a fixed value to 
symmetry (zero gradient) at X1=TIMESW for axisymmetric 
fluted mixers 

0 Default, suppress printout of input cards 

1 Print input cards 

0 Default, automatically computes number of columns (for 

computation of U2 from continuity equation solver) from 
rectangular grid input 


n For axi-flows (NAXI=1), assumes n radial lines for computation 
of secondary velocities 


NAXI (V) 
NBC (I) 


NELEM (IV) 


0 Default 

1 For axi-flows 

0 Default 

n Numbers of nodes for which gradient boundary condition 

(KBNOTW=l ) will be used 

0 Default, automatically computes number of elements assuming 

rectangular grid 

1 For axi-flows (NAXI=1), node numbers for each element will 

be input, see Section 4 
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NELEMS (IV) 

0 

Default 


n 

Number of elements to delete when boundary condition changes 
at X1=TIMESW for axisymmetric fluted mixers 

NEQKNN (I, II, III, IV) 

1 

Default 


n 

Number of dependent variables to be integrated in XI direction 

NNODE (IV) 

0 

Default, number and location of nodes automatically computed 
assuming rectangular, grid 


1 

For axi-flows (NAX1=1), node locations will be individually 
input, see Section 4 

NNODES (IV) 

0 

Default 


n 

Number of nodes to delete when boundary condition changes 

NPVSX (I,II,III,IV) 

n 

Number of entries in pressure table 

NSCX (I,II,III) 

0 

Uniform X3 interval in rectangular discretization (Default) 


1 

Nonuniform X3 interval in rectangular discretization 

NSCY (I,1I,III) 

0 

Default, uniform X2 interval in rectangular discretization 


1 

Nonuniform X2 interval in rectangular discretization 

NE1E2 (I) 

0 

Default, laminar flow 


1 

Turbulent flow (see Chapter IV for input for various turbulence 
models) 

NU2POS (II) 

n 

Number of entries in table for variable geometry in X2 
direction 

NU3POS (II) 

n 

Number of entries in table for variable geometry in X3 direction 
(Note: even if variable geometry is only being used in one 
coordinate, NU2POS and NU3POS must both be input as the 
same value) 

This part of Section 3 for inputting integer variables is ended by a card with SEND starting in Column 
2. 

The second part of Section 3 consists of inputting noninteger input variables. It begins with the 
card SNAME02 in Column 2, followed by the input variables in the same format as NAMEOl , and 
ended by the card SEND starting in Column 2. A complete list of all the possible input variables is 
given in Table XII, which was computed for Case I. For the four test cases presented here, the required 
variables and their input units are; 

AO (III) 


Area of primary flow (F^) 

A1 (III) 


Area of secondary flow (F^ ) 

A2 (III) 


Area of diffuser entrance (F^ ) 

A3 (III) 


Area of diffuser exit (F^ ) 
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ADUCT (IV) 


ALC-(IV) 


BMDOT (III) 
CMDOT (III) 

DELP (I, II, III, IV) 
PCNT (I) 

PMDOT (III) 

REEL (I,II,III,IV) 
RHOINF (III) 

SMDOT (III) 

TD (I,II,III,IV) 

TO (I,II,III,IV) 
TIMESW (IV) 

TOFINF (I,II,III,IV) 
UB (III) 

UC (III) 

UINF (I,II,III,IV) 

UP (III) 

VSTART (I,II,III,IV) 

X2SCLE (II,III) 
X3SCLE (II,III) 
X2SHFT (II) 


Multiplier for confined flow pressure computation. Total area 
of cross section = ADUCT times the computed area of the 
input cross section as determined from the sum of all the 
input elements. ADUCT has a default value of 2.0 (Cases 
II and III). 

Program computed length scale for non-dimensionahzation. 
For rectangular grids, program automatically computes as 
shortest distance between adjacent X2 or X3 values of the 
nodes. For axi-flows (NAXI=1), user input as shortest side 
of the finite element triangles. 

Flow rate of boundary layer control nozzle (Ibm/s) 

Flow rate of Coanda nozzle (Ibm/s) 

Percent increment of TD at which output is desired 

Percent difference from edge velocity for computation of 
boundary layer thickness 

Flow rate of primary nozzle (Ibm/s) 

Reference length (F) 

Reference density (Ibm/F^). If not input, computed from 
reference temperature and first entry in pressure table. 

Flow rate of secondary flow (Ibm/s) 

Length of XI solution, starting at TO(F) 

Initial XI station (F) 

XI location (in feet) where boundary condition changes for 
axisymmetric fluted mixers. Default = TO+2*TD 

Reference stagnation temperature (°K) 

Velocity of boundary layer control nozzle (F/s) 

Velocity of Coanda nozzle (F/s) 

Reference velocity (F/s) 

Velocity of primary flow (F/s) 

Percent increment of TD at which transverse velocity (U2) 
computation starts - Default = 1 00% 

Multiplier on X2 to convert rectangular discretization to feet 

Multiplier on X3 to convert rectangular discretization to feet 

Shifts X2 coordinates by an amount X2SHFT(F)/ALC, where 
ALC is the program computed length scale (see entry 3 
in Table XI 
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After the END card for NAME02, Section 3 is ended by the card FEDIMN starting in Column 1 
calling the subroutining FEDIMN which computes dimensionalizing quantities, vector lengths, etc. 
based upon the variables input in NAMEOl and NAME02. 

Section 4. 

This section computes the node location and finite elements. It begins with the card LINKl in 
Column 1 and 1 in Column 10. Three modes may then be used. 

Mode 1 : Automatic Uniform Rectangular Discretization 

Occurs for NSCX = NSCY = 0 (Case II), where user set X2SCLE = desired element width in the X2 
direction and X3SCLE = desired element height in the X3 direction. On the next card in any column 
(i.e., input for this card is in free format), indicate the number of the first node in the X2 direction 
(always 1 ). Skip one or more places and add the last number in the X2 direction. Place a comma 
immediately after the last X2 node number. Repeat for X3 direction, e.g., if we desired 2 columns 
of 1 1 nodes in the X2 direction, the card would appear as: 

1 11 , 1 2 , 

The sequence is completed by placing a “T” either on this card or a following card. 

Mode 2: Automatic Nonuniform Rectangular Discretization 


Occurs for NSCX=NSCY=1, cases (I, III). Set X3 discretization first; X2 second. Data are input in 
sets of 3 integers separated by a comma after the third integer. The first integer identifies the interval 
number between nodes, the next two indicate interval width (or height) as a dimensionless ratio. The 
first integer in successive groups need not be consecutive, i.e., if the first 5 intervals are to be the same 
distance apart, then it is allowable to skip from ! to 5 for the first integer. For example (see Table I), 
there is one interval of width 1/60 in the X3 direction and 20 intervals of height 4/1200, 4/1200, 
4/1200, 8/1200, 1/120, 1/120, 1/60, 1/60, 1/60, 1/60, 1/30, 1/30, 1/30, 1/30, 1/30, 1/30, 1/30, 

1/30, 1/30, 1/30, 1/30 in the X2 direction. Hence, the grid is 21 nodes high (X2) by 2 nodes wide 
(X3). These dimensionless intervals are converted to feet by being multiplied by X2SCLE and 
X3SCLE. 

Mode 3: Nonuniform, Nonrectangular Discritization 


Occurs for NNODE=l, NELEM=1, in NAMEOl (Case IV). For this mode, the three node numbers 
for each element are entered in integer form (free format with no commas) ended by a “T”. On the 
next card, in Column 5, put a 1 . On successive cards (in Format 2E20.8), input the values of X3 and 
X2 for each node, starting with node number 1 . End this input with any large number greater than 
10^°. The numbering system used to generate the finite elements and node locations for Case IV is 
shown in Figure 32. 

Section 5. 

Starting in Column 1 COMTITLE, which designates that the following card will be a title card 
(Columns 1-80) printed on the cover page of the output. End this section with a DONE card, i.e., 
starting in Column 1 , the letters DONE. 
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Section 6. 


Starting in Column 1 DESCRIPT, then in Columns 11-13 the integers 204 followed by one or more 
spaces and a “T”. Note that any message may be placed on this card after the “T” without affecting 
the data. Up to ten successive title cards may then be added which will be printed at the head of each 
output print station. End this section with a DONE card. 

Section 7. 

Starting in Column 1 DESCRIPT, then in Columns 1 1-1 3, the integers 332 followed by one or 
more blanks and a “T”. The following cards are the output heading descriptors which are printed at 
the beginning of each print station. This section is ended by a DONE card. 


n _ _ x! o 

dCCllUll O. 

Starting in Column 1 lONUMB, followed by -1 in Columns 9 and 10. The next cards print the 
appropriate variables from NAME02 (see Table XI) in the spaces provided (5 per output line in free 
format) in the output heading descriptors of Section 7: A “999” gives a line of x-x’s across the page, 
while a “200” leaves the space blank. This section is ended by a “T”. 

Section 9. 

Starting in Column 1 MPARA, followed by -1 in Columns 9 and 10. These cards multiply each 
input variable from Section 8 by the corresponding values (mostly conversion factors) from NAME02. 
Note that entry 2 in NAME02 multiplies the variable by 1 .0. This section is ended by a “T”. 

Section 10. 

Starting in Column 1 DESCRIPT, then in Column 1 1 the number 203, followed by one of more 
blanks and a “T”. The succeeding cards give the titles to be printed above the output variables (inte- 
grated or computed). The format is 5H16 (5 variable titles per card). This section is ended by a 
DONE card. 

Section 11. 

This section designates which variables are to be printed out at all node point locations. Any vari- 
able ending with the numerals 248 represents a dependent variable, those with a 249 represent the 
derivative of the dependent variable. The numberals preceding the 248 or 249 designate which de- 
pendent variable is to be printed. See Section 16, for a list of the corresponding integer value for each 
dependent variable. The following computed variables may also be printed out at each node. 


Location 

Variable 

Location 

Variable 

285 

Static temperature 

1247 

Effective viscosity 

320 

Static enthalpy 

334 

Effective Prandtl number 

284 

Density 

314 

Effective Schmidt number 

278 

Cp 

292 

Laminar viscosity 


This section is ended by a “T”. 
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Section 12. 


Starting in Column 1 lOMULT, followed by -1 in Column 8. The following cards multiply the 
input variables from section 1 1 by the appropriate values in NAME02. For division, put a minus 
sign in front of the integer value. This section is completed with a “T”. 

Section 13. 

Starting in Column 1 COMOC. This card prints the COMOC logo at the beginning of the output 
print. 

Section 14. 


Starting in Column 1 DESCRIPT. The following title cards (Columns 1-80) are printed once at the 
start of the run. A maximum of four title cards is allowed. This section is completed with a DONE 
card. 

Section 15. 

Starting in Column 1 VX3ST designating that the next card will give the XI locations (in feet) at 
which the pressure will be input. If the 3DBR option is used (Case I), this pressure profile will be 
used for the complete calculation. If the 3DPNS option is used (II, III, IV), this pressure profile is 
used to specify the initial pressure and to compute the initial pressure gradient. The XI locations are 
input in free format separated by one or more spaces and are terminated by a “T”. The next card has, 
starting in Column 1 VPVSX, followed by a card with the corresponding pressure values (in PSFA) 
terminated by a “T”. 


Section 16. 


Starting in Column 1 IPINT and in Column 1 1 put -1 . The following card designates the integer 
array numbers of the dependent variables. The first NEQKNN of these values are integrated and also 
U2. Starting in any column, list the integer values of the dependent variables which correspond to 
the following key: 

Integer Variable 


1 

2 

3 

4 

5 

6 
7 


U1 

U2 

U3 

H 

k 

d 

Open 


8 

9 

10 


yHa 

yN2 


A free format is used, the sequence being terminated with a “T”. 
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Section 17. 


Boundary conditions are required at each boundary node point for all integrated, dependent vari- 
ables. The following input instructions refer to any dependent variable. There are three modes under 
which the boundary condition may be defined: (1) zero normal gradient, (2) value fixed to its initial 
value, (3) normal gradient defined as a finite value, i.e., the gradient boundary conditions. 

Mode 1 : Zero Gradient Boundary Condition 

This represents the default boundary condition, hence any boundary node whose boundary condi- 
tion for a particular dependent variable is not input will automatically be given a zero normal gradient. 

Mode 2: Boundary Node Fixed to its Initial Value 


This mode may be input in two ways. First, for rectangular grids, we define the boundaries of the 
cross section as TOP, BOTTOM, LEFT, and RIGHT. The code has been programmed so that the fixed 
boundary condition for the nodes along a given boundary are not required to have the same boundary 
condition. To use this option (see Cases I, II, III) on the first card, put KBNO in Column 1 , then in 
Column 1 1 put the integer number of the dependent variable whose boundary condition is to be de- 
fined. The next card may have the general image as follows: 


Column = 1 

1 1 2 1. 

31 

41 

5 1 

61 

7 1 

Card = BOTTOM 

NB TOP 

NT 

RIGHT 

NR 

LEFT 

NL 


This general card designates that the first NB nodes of the bottom row, NT nodes of the top row, NR 
nodes of the right column, and NL' nodes of the left column will be held constant and equal to their 
input values. By leaving out any of the integers (NB, NT, NR, NL) all the nodes on that particular 
boundary are held constant. For those nodes not defined, the default boundary condition (Mode 1) 
applies. This input is terminated by a DONE card. 

The second option for this mode is to pick specific nodes in the cross section to be kept fixed to 
their input values (as in Case IV), or to be recomputed by some algebraic relation (as in Case I). The 
input for this option is on the first card KBNO in Column 1 followed by the integer number of the de- 
pendent variable in Column 1 1 . The next card is ADD starting in Column 1 . The following card has 
the integer numbers of the nodes in free format separated by blanks and terminated by a “T”. This 
option is then completed by a DONE card. 

Mode 3: Gradient Boundary Condition 

This option is used when KBNOTW=l . It sets the derivative of U1 at the wall equal to the wall 
shear stress (computed in WLFLXS) divided by the laminar viscosity at the wall. The input form (see 
Case I) is as follows. On the first card, starting in Column 1 , IBORD. The next card, again in Column 
1, BOTTOM, -TOP in Column 21 , and DONE in Column 41 . These two cards orient the elements so 
that the gradient boundary condition may be properly applied. For the next card, starting in Column 
1, KBNO, followed on the same card by a 1 in Columns 10 and 12. For the final card, starting in 
Column 1, BOTTOM, 0 in Columns 11, 13, and 15, 0.0 in Column 18, 2 in Column 23, 0.0 in Column 
27, 2 in Column 34, 0 in Column 39, and a “T” in Column 41. 
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Section 18. 


This section contains the format for inputting the initial values of the variables at each node point. 
The format for each variable is as follows: 

Stagnation Temperature (Case III) 

On the first card, starting in Column 1 , VTEMP. If the stagnation temperature is input in °R, then 
in Column 1 1 put a -58. This divides the input values by the 58th entry in NAME02 which is 
TOFINF, the reference temperature in °R. The next card has the values of the stagnation temperature 
for each node. The order that the values appear on the card(s) correspond to the order of the node 
numbers. For Modes 1 and 2 of Section 4 (Rectangular Grid), this order proceeds from the bottom 
node of the left column across the bottom row, then to the next row (left to right), etc. For Mode 3 
of Section 4, the node numbering follows from the order of the cards for which the values of X2 and 
X3 are prescribed. Repeated values may be defined as n-A, where the n successive nodes would each 
have a value of A. This data set is terminated by a “T”. If this data set is not prescribed, then each 
node is assumed to have an input stagnation temperature of TOFINF. 

Integrated Variables (Including U2) 

The input for these variables is as follows. Starting in Column 1 VYY, then in Column 1 1 the 
appropriate nondimensionalizing factor (-27 for velocities, etc.). Note, that if Columns 1 1-14 are left 
blank, then the input values are nondimensional. The succeeding cards contain the values of the vari- 
able at each node in the same format as the stagnation temperature, terminated by a “T”. The final 
card has a VYYEND in Column 1 , followed by the integer value n in Column 1 1 , where n denotes the 
appropriate integer (see key in Section 16) for each variable. 

Section 19. 

This section sets the laminar viscosity equal to a constant value at each node point, see Case II. 
Starting in Column 1 put VSUTHLD, on the following card in any column after 1 put 4*0.0 then one 
or more blanks, and finally the constant value in free format. This section is terminated by a “T”. 

Section 20. 

This section inputs a turbulent viscosity at each node point, see Case III. On the first card, starting 
in Column 1, VEPSILON, followed by -38 in Columns 11-13, if the turbulent viscosity is to be input 
in dimensional units (LBM/F-S). If input is nondimensional form, leave Columns 1 1-13 blank. On the 
next card input the values for each node, in the same order as for the velocities, etc. Terminate this 
section with a “T”. 

Section 21. 

This section inputs the table for a variable cross-sectional area solution domain in the X2 direction 
(see Summary of Key Items for Input for a more general description of the input options for variable 
area). On the first card, starting in Column 1 VU2POS. The next card contains the XI locations, in 
feet, where the values of X2 are to be input. The number of input values is equal to NU2POS in 
NAMEOl . They are input in free format, and terminated with a “T”. On the next card, starting in 
Column 1 , VU2VAL. The following card contains the values of X2, in feet, at the top of the solution 
domain for each corresponding XI station. The final card contains the values of X2 at the bottom of 
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the solution domain, terminated by a “T”. Variable geometry in the X3 direction is similarly input 
using VU3POS and VU3VAL cards. 

For axisymmetric fluted mixer flows (NAXI=1), the variable geometry in the radial direction is in- 
put as follows. The first card is VU3POS, followed by a card 'with the XI locations (in feet) for 
specifying the inner and outer radii of the duct, the number of entries being equal to NU3POS in 
NAMEOl. They are entered in free format and terminated by a “T”. The next card contains VU3VAL 
starting in Column 1. The following two cards contain the inner and outer radii, respectively, of the 
duct (in feet) at the corresponding XI locations. The last card is terminated with a “T”. 

Section 22. 


This section sets up the node locations for computing the U2 and U3 velocities along radial lines 
from the continuity equation in polar coordinates (N.AJKI=1, see Case IV). On the first card, starting 
in Column 1, CNTPTS, followed by -1 in Column 1 1. The next card has the integer number of nodes 
in each radial line (there are LG entries) completed by a “T”. The next card, starting in Column 1, is 
CNTNDS, followed by -1 in Column 1 1. The final cards contain the integer numbers of the nodes 
along each line, terminated by a “T”. See Figure 32 for the notation used to obtain the numbering 
for Case IV. 


Section 23. 


For the first card, starting in Column 1, LINKCALL, followed by -1 in Column 1 1 . The next card 
describes the calls to various subroutines performed at each integration step. The calls are inputted as 
pairs of two integers followed by a comma. The end of the calling sequence is determined by the 
single integer 0, followed by one or more blanks and a “T”. The purpose of the subroutine called by 
each integer pair is as follows: 


Integer 

Pair 

Subroutine Called 

Purpose 

2 

4 

CONTES 

Computes U2 Velocity from continuity equation 
knowing U1 and U3 (NAX1=0). Or, computes 
U2 and U3 from continuity equation in polar 
coordinates knowing U1(NAXI=1). 

5 

6 

SETDIF 

Computes diffusion coefficients 

2 

15 

TBRTHK 

Computes integral parameters for boundary layer flows 

2 

9 

DRHOBL (IGAS=0) 
DRH0GSAIGAS=1) 

Computes nodal values of density 

2 

3 

WLFLXS 

Computes skin friction and wall shear stress 

1 

5 

DPDX 

Computes pressure gradient, as either input 
(3DBR) or calculated (3DPNS) 

1 

6 

SOURCE 

Computes source term for each dependent variable 

2 

19 

WLSTAR 

Computes k&d at 1st node off the wall 

5 

1 

NWGEOM 

Computes new coordinates for variable geometry 

5 

11 

OUTPAR 

Computes ejector performance parameters 
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The next card has starting in Column 1 LINKS, followed by a 4 in Column 10. This calls the sub- 
routine DIMEN which finalizes all the nondimensionalization. The final card has starting in Column 
1 , LINKl , followed by a 3 in Column 1 0. This card establishes which nodes for each variable are to be 
integrated as well as setting up the coordinates in the proper nondimensional form. 

Section 24. 

Starting in Column 1 , QKNINT. This card directs the program to the integration loop. 

Section 25. 

Starting in Column 1 , DESCRIPT, followed by a 3 in Column 10. The following cards denote the 
printout which will appear at the end of the program output. This section is ended with a DONE card. 

Section 26. 

Starting in Column 1 , COMOC. This card prints the COMOC logo after the output of Section 25. 

Section 27. 

Starting in Column 1 , EXIT. This tells the program that the output is complete and to end execu- 
tion. 


OUTPUT DESCRIPTION 

The previous sections dealt with preparing the input data decks for Cases I-IV shown in Tables. 

Ill - VI. These data decks produce the output format at each print station as shown in Tables VII - X. 
The output resulting from pertinent sections of the input data decks has been labeled on these tables. 

SUMMARY OF KEY ITEMS FOR INPUT 

This section will summarize and provide a more complete description for inputting some of the 
more important aspects of the flow systems under investigation. 

1 . Turbulence Modeling 


The COMOC code currently has five different turbulence models available to the user. A 
description of each as well as input instructions follows. 

1-a. Mixing Length Theory: This model computes the turbulent viscosity from classical mixing 
length theory, including von Driest damping in the near wall region and intermittency in the outer 
freestream region. This option is used by setting NEIE2=1 , in NAMEOl . For this model k&d are not 
to be integrated and hence should not be included in IPINT (Section 16) or in the determination of 
the number of variables to be integrated (NEQKNN). 

1-b. Two-Equation (k&d) Model with Initial Profiles for k&d Computed from Mixing Length Theory 
This option is used in case I. Here the turbulent viscosity is obtained at each downstream station by 
integrating the k&d equations. However, the initial profiles are obtained by using mixing length theory 
at the initial XI station to compute the mixing length and turbulent viscosity. Algebraic relations are 
then used to compute the initial k&d profiles. To use this option input NE1E2=1 in NAMEOl , include 
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k&d when determining the number of equations to be integrated (NEQKNN) and finally by including 
5 and 6 in the first NEQKNN entries of IPINT (Section 16). 

1-c. Two-Equation (k&d Model with Initial Profiles for (k&d) Input at Each Node Point: For this 
option the initial profiles for k&d are input following the format in Section 1 8. If k&d are to be input 
in dimensional form (F^/S^ and F^/S^) then the nondimensional factors for Column 1 1 on the VYY 
cards' are, -90 and -89, respectively. Also required are that NE1E2=1 and ITKE= 1 in NAMEOl . 
NEQKNN and IPINT (Section 1 6) include k&d as for Case 1-b. 

1-d. One-Equation (k) Turbulence Model: For this case the turbulent viscosity is obtained from k 
(integrated) and a mixing length C computed algebraically for free shear flows. For this case the initial 
k profile is input as in Case 1-c above. In NAMEOl include NE1E2=1, ITKE=1, and IFSL=1. Also 
include k when determining NEQKNN and hence in IPINT. 

1-e. Constant Turbulent Viscosity Model: This option was used in Case III. Here the turbulent 
viscosity at each node is input in VEPSILON (Section 20). No additional input in NAMEOl is re- 
quired and k&d should not be included in NEQKNN or IPINT. 

2. Boundary Conditions for Turbulent Flows 


Turbulent flows are characterized by steep gradients in the near wall region. For this reason 
it was found to be advantageous to replace the no-slip wall velocity boundary condition for Ul, and 
to calculate the values of k&d at the node off the wall from algebraic relations. This requires fewer 
node points in the near wall region and hence allows faster running times. The gradient boundary con- 
dition is used by inputting KBNOTW=l, and NBC=n where n is the number of columns, in NAMEQl, 
and by using Mode 3 of Section 17 for the Ul boundary condition at the wall. The computation of 
k&d off the wall is done algebraically by calling 2 19 in LINKCALL (Section 23). Hence, the nodes 
off the wall for k&d need not be integrated and may be eliminated (Mode 2 - second option, of 
Section 17). 

1 

3. Variable Geometry 


Rectangular or axisymmetric flows whose cross section does not remain constant in the XI 
direction may be computed using a coordinate transformation. For rectangular flows, the required 
input for this option is as follows. In NAMEOl , set NU2POS=NU3POS==n, where n is the number of 
entries in the table of XI vs X2 and/or X3. In the LINKCALL card (Section 23) include the integer 
pair 5 1 , which calls the subroutine for recomputing the coordinates at each XI station as well as com- 
puting the additional terms in the differential equation resulting from the coordinate transformation. 
For rectangular flows, two options are possible. The first treats variable geometry flows which are 
symmetrical about the centerline. For this case in NAME02 put ADUCT=2.0 (default). The flow is 
then calculated in the bottom half of the cross-sectional solution domain only. The input cards with 
the XI vs X2 and/or X3 values are read in as in Section 21 (see Case II). For nonsymmetrical flows, 
the entire cross-sectional domain must be included. This option is set by putting ADUCT=1 .0 in 
NAME02. The XI vs X2 and/or X3 tables are then read in as before. Variable geometry for axi flows 
whose inner and/or outer radii vary with XI may likewise be calculated using the input described in 
Section 2 1 . 


I 
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4. Ejector Parameter Output 

For ejector/diffuser systems, several parameters are of interest, i.e., thrust augmentation, skew- 
ness, etc. (see Appendix B for the definition of these terms). These parameters are computed and 
printed out at each output station as follows. The enlarged output print of Section 7, 8, and 9 is used 
(see Case III). In the LINKCALL (Section 23) include the integer pair 5 11, which calls the sub- 
routine to compute the parameters. Also requireid in input are AO, Al, A2, A3, UB, UC, UP, BMDOT, 
CMDOT, PMDOT, and SMDOT in NAME02. 

5. Node and Element Input for Axisymmetric Flows 


For axisymmetric flows. Mode 3 of Section 4 is used. This requires NAX1=1 , NNODE=l , 
NELEM^l in NAMEOl, and the specification of ALC in NAME02. Also, for the computation of U2 
and U3 from the continuity equation in radial coordinates. Section 22 giving the node numbers along 
radial lines is required (see Case IV). 

6. Changing Boundary Condition for Axisymmetric Flows 

At the downstream X location where the inner wall (conical center) ends, the boundary con- 
dition for those nodes that were along the wall changes from no slip to symmetry, i.e., zero gradient. 
Also at this station, those finite elements and nodes that were along the wall require further exam- 
ination. Figure 33 shows these elements and nodes at the initial X station and at the station where 
the inner wall ends. The nodes are depicted by the numbers 31 to 40, and the elements by the circled 
numbers from 49 to 56. As can be seen, at the X station where the inner wall ends, the elements 
50, 52, 54 and 56 shrink to zero area. This causes computational problems and therefore these ele- 
ments must be deleted from the computational domain at this X station. Also, nodes 36, 37, 38, 39, 
and 40 become physically identical. This will not cause problems as each node is computed indepen- 
dently. However, since element 56 has been deleted, node 40 is no longer used and must also be 
deleted. 

The input instruction to perform these procedures are as follows: 

6-a. The variable geometry for axisymmetric flows (NAX1=1 ) is calculated by inputting the inner 
and outer radii as functions of XI as per Section 21 and by calling the subroutine NWGEOM by in- 
putting the integer pair 5 1 in the LINKCALL as per Section 23. 

6-b. The Xl station in feet where the inner wall ends and hence the boundary condition changes is 
input as the value of TIMESW in NAME02. 

6-c. The number of initial nodes whose boundary condition changes at X1=TIMESW is input as 
the integer KBNOS in NAMEOl. For the example shown in Figure 33, KBNOS=5. It should be noted 
that the last KBNOS nodes change their boundary condition input as per Section 4 (Mode 3). There- 
fore, care must be taken when setting up the node locations to insure that these nodes are read in last. 
Also, even thougli node 40 is deleted at Xl^TlMESW, it is still included in the determination of 
KBNOS. 

6-d. The number of elements to be deleted at X1=T1MESW is input as the integer NELEMS in 
NAMEOL For the example in Figure 2, NELEMS=4. Again, the last NELEMS elements input as per 
Section 4 (Mode 3) are deleted. 
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6-e. The number of nodes to be deleted at X1=TIMESW is input as the integer NNODES in 
NAMEOl (NNODES=l for the example in Figure 33). Again, the last NNODES nodes of the input 
are deleted. 


CONCLUDING REMARKS 

This manual is intended to provide the user who is not familiar with COMOC with a step-by-step 
outhne of how to prepare an input data deck. The test cases presented were chosen for two reasons. 
First, they represent typical problems when dealing with ejector/diffuser or fluted nozzle mixer sys- 
tems. Secondly, they contained all the pertinent options currently in effect in the COMOC code. 
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APPENDIX B 


DEFINITION OF EJECTOR/DIFFUSER PARAMETERS 
1 . Ejector Inlet Area Ratio 

The following discussion presents definitions of ejector terminology with respect to ejector 
station locations shown in Figure B-1, which along with the references used herein may be found at 
the end of Appendix B. 

There is no ambiguity in the definition of inlet area ratio, Aj/Aq, when the ejector consists of a 
single primary nozzle discharging axially into a symmetric mixing chamber, and the primary nozzle 
is assumed to be ideal (von Karman, 1969; Bevilaqua, 1973). Under these conditions, the inlet area 
ratio is defined as the ratio of the secondary stream cross-sectional area at the nozzle exit plane, to 
the primary nozzle exit area. 

The first difficulty arises when the performance of an actual primary nozzle is considered. If 
the primary flow is subsonic and the nozzle is conical, the primary stream will initially contract 
until there is a static pressure balance between the primary and secondary flows. The use of inlet 
area ratio as a performance correlation parameter arises from a one - dimensional analysis in which 
the static pressures of the two streams are assumed to be equal over the areas Aq and Aj (e.g., Camp- 
bell & O’Hain, 1967). Better correlation could be expected, therefore, if the primary nozzle effective 
area were used rather than the geometric area (Salter 1972). McClintock and Hall Hood, (1966), in 
discussing the performance of ejectors with supersonic primary flows, suggested that the fully ex- 
panded area should be used for ejector performance calculation. It is, however, more customary to 
define the inlet area ratio as a geometric area ratio, and to include supplementary curves of primary 
nozzle performance (Campbell, et. al, 1973). 

When the primary flow consists of multiple jets discharging into the ejector at a number of dif- 
ferent locations, the task of correlating performance on the basis of flow areas at some station where 
the pressure is assumed to be equalized would be a formidable task. It would be necessary to deter- 
mine the expansion of each jet to the common reference plane, and to calculate each jet velocity at 
that plane. Under these conditions Quinn (1972) chose to define inlet area ratio as 


A\ _ Total duct area at the primary exit plane 

Aq Sum of geometric primary and boundary layer nozzle areas 


From the practical aspects of ejector design and performance evaluation, inlet area ratio should be 
defined in terms of geometric areas, and Quinn’s definition for the case of multiple injection appears 
to be the most appropriate. We suggest, therefore, that inlet area ratio should be defined as 


lAR = 


Ejector duct area at primary nozzle exit plane 

Total geometric area of primary, Coanda and Boundary layer nozzles 


2. Diffuser Area Ratio 


Ambiguity could arise here too if an attempt were made to define the ratio in terms of areas 
normal to the flow directions. For example, at the diffuser exit under conditions where the diffuser 
doors are deflected backward, one might define the exit area as the plane normal to the diffuser 
centerline through the trailing edge of the forward door (Bell, 1973). This definition results in the 
diffuser area ratio being a function of the deflection angle even when the diffusion angle is kept con- 
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stant. The reason for this definition was an attempt to correlate ejector performance independent of 
deflection angle. For the purposes of design and performance stimation, it is more desirable to ex- 
press the diffuser area ratio in terms of the physical diffuser inlet and exit boundaries. We suggest 
that diffuser area ratio should be defined as: 


DAR = 


Cross-sectional area at Station 3 
Cross-sectional area at Station 2 


(3) 


3. Gross Thrust 

Gross thrust is the thrust produced by the flow leaving the ejector. A simple integration of the 
momentum flux over the diffuser exit plane area (Station 3) would not yield the correct gross thrust 
value because the stream lines may be diverging and the pressure may be varying over this area. The 
most direct way to determine the gross thrust value is to integrate the momentum flux over the sur- 
face where the static pressure is everywhere equal to the ambient pressure, i.e., over Station 3'. 


If n3' is the unit vector normal to the surface 3' then 


Gross Thrust == / 


A3' ^3' (P3'U3') * n3' da 


(4) 


4. Momentum Drag 


Momentum drag is the loss in momentum of the freestream due to the presence of the ejector. 

Momentum drag = fniVoo = Pl^i * nj da (5) 

5. Net Thrust 

In the case of a propulsion unit whose exhaust is aligned such that the gross thrust acts along 
the fli^t path, the net thrust is defined as the difference between the gross thrust and momentum 
drag. When the gross thrust acts at an angle to the freestream direction, the resultant force may be 
computed by vector subtraction. It is suggested that, rather than computing net thrust, the magni- 
tudes and directions of the gross thrust and ram drag should be given. Pitching moments may then 
be computed from the sum of the two moments about any desired point. 

6. Mass Flow and Mass Flow Ratio 

The mass flow rate through any surface is obtained by the integration of the components of 
(pU) normal to the surface, over the surface area. 


m = /^(pU) • nda 


( 6 ) 


The total primary mass flow rate is (irip + m^). The question is now raised as to whether or not the 
boundary layer control flow should be included in the definition of mass flow ratio. From a practical 
point of view, all air supplied to the ejector jets must be supplied by power plants on board the air- 
craft. If this air were not used in an ejector, it could be employed to produce useful thrust. For 
example, if the boundary layer control air were extracted from the engine bleed ports, there would 
be a reduction in engine performance, either in the form of a direct thrust loss, or in an increase in 
fuel consumption if the power lever were advanced to maintain the same thmst. It is reasonable, 
then to consider all of the injected mass flow in the definition of mass flow ratio. In addition, this 
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definition would be consistent with the. inlet area ratio definition. From these considerations, it is 
recommended that the following definition should be used: 


Mass Flow Ratio 


Total ejector discharge mass flow rate 
Total injected mass flow rate 


7. Thrust Augmentation Ratio 


(7) 


Von Karman (1949) defined thrust augmentation ratio as the ratio of the actual thrust pro- 
duced to the primary jet momentum with the nozzle installed in the ejector. He commented that 
this is “somewhat arbitrary.” Campbell and O’Hain (1967) and all subsequent ART literature, de- 
fine thrust augmentation ratio as the ratio of ejector thrust to the isentropic primary thrust that 
would result if the same primary mass flow rate were discharged to atmospheric pressure, from the 
same stagnation conditions. Since the local pressure at the nozzle exit plane is lower than ambient 
pressure, this definition implies a larger nozzle area when discharging to ambient pressure. The 
reasons for this definition are two-fold. First, a certain mass flow rate at known stagnation condi- 
tions is provided by the aircraft power plant. If the engine is to run at its design point when opera- 
ting in conjunction with an ejector, the primary nozzles must be sized correctly to allow the same 
mass flow through the system. ARL’s definition is, therefore, the ratio of ejector thrust to that 
which could be produced by the same engine exhausting through an isentropic nozzle. The second 
reason for adopting this definition is its convenience. The ejector thrust is measured during a test 
and the isentropic thrust is computed from the stagnation conditions and mass flow. Another 
definition used by Boeing/NASA (Campbell 1973) is the ratio of ejector thrust to actual primary 
thrust with the ejector shrouds removed. In this approach, the primary thrust is reduced due to the 
lower mass flow (unless the nozzles are choked), and the primary nozzle thrust efficiency is included. 
This definition leads to a larger value of thrust augmentation ratio than ARL’s. Alperin and 
Marlotte (1971) discuss various definitions and conclude that a parameter based on engine power is 
the most appropriate. Their definition consists of the ratio of ejector thrust to reference engine 
thrust based on the same reference engine power. As indicated by Alperin, a problem with this 
definition is that it requires an assumption as to the reference engine nozzle area. 


It is our recommendation that ARL’s definition of thrust augmentation should be used, and 
that all jets dischargina into the ejector should be considered in the evaluation of the denominator. 
Denoting Vj^ as the isentropic velocity of a jet, which can be calculated from its stagnation pressure, 
stagnation temperature, and ambient pressure, then thrust augmentation ratio 0, is defined as: 

j _ Ejector Gross Thrust 

2 (AVi,)i i = P.=.b (8) 

No attempt will be made to define net thrust augmentation ratio, because of the lack of a suit- 
bale definition for net thrust. 


8. Propulsion Efficiency 

The basic definition of propulsive efficiency is the ratio of output work rate to input energy rate. 
The output work rate consists of that due to the component of gross thrust in the direction of flight 
and that due to momentum drag. Thus, 

Output work rate = Fq '^oo (9) 

— y — ^ 

where F q is the gross thrust vector and is the momentum drag. 
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Input energy rate is the sum of the input energy rates of all of the jets. To be consistent with 
our definition of thrust augmentation ratio, it is recommended that isentropic jet velocities should 
be used for the calculation of input energy rates. Then propulsive efficiency is: 


’^PR 




1/2 Z (mVfgjj 


i = c, p, b 


( 10 ) 


9. Energy Transfer Efficiency 

This is the ratio of the rate of energy addition to the freestream by the ejector to the input 
energy rate of the jets, i.e.,: 


^TR " 


m'j’|V3' 





i = c, p, b 


where V3', is the integrated average velocity in Station 3. 


( 11 ) 


10. Skewness 

Skewness is a convenient means of expressing the “flatness” of the velocity profile at any station. 
It has been denoted in the literature by X (Von Karman, 1949), Kj^ (London, 1955 and SAE, 1960) 
and by (5 (Quinn, 1972, and Bevilaqua, 1973). 

Skewness is defined as the ratio of the actual integrated momentum to the bulk average momen- 
tum at the station under consideration. Von Karman was the first investigator to demonstrate the 
importance of this parameter on ejector performance. He showed that high inlet skewness was bene- 
ficial, and suggested that the high efficiency of the Coanda ejector was largely explained by this 
phenomenon and “has little to do with the mysterious flow of the fluid around sharp edges.” 

Bevilaqua (1974) has demonstrated the importance of low values of skewness at entrance to the 
diffuser. 


The bulk average momentum at any station is 
mV = ^ ■ n da 

— ^ ^ 

and the actual momentum is (pU) • Uda. Skewness is defined, therefore, as 

A /a (pU) • Uda 
J A U ’ n da 

1 1 . Loss Factor 


( 12 ) 


(13) 


This is defined for any component as the ratio of total pressure loss across the component to the 
dynamic pressure at entrance to that component. Thus, s 



( 14 ) 
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where AP is the loss in total pressure, and q is defined by 

q = -^ (pU) • Udo at entry 

1 2. Ideal Isentropic Nozzle Thrust 

This parameter was discussed under “Thrust Augmentation Ratio,” where the following defini- 
tion was recommended: 

Fis - S (mVis)i i = p,c,b (15) 

i 

1 3. Actual Isentropic Nozzle Thrust 

This is defined as the sum of the isentropic thrusts of all of the jets discharging into the ejector. 
One of the program output parameters will be a complete description of the pressure field in the 
ejector. The local pressure at each nozzle exit location will be used, in conjunction with the stagna- 
tion pressures, temperatures, and mass flow rates, to calculate the actual isentropic nozzle thrusts. 
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Table [ 

Ejector Results at x = 0.0, 15.2, and 27.9 cm (cont) 
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Table I 

Ejector Results at x = 0.0, 1 5.2, and 27.9 cm (cont) 
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Table I 

Ejector Results at x = 0.0, 15.2, and 27.9 cm (cont) 
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Table I 

Ejector Results at x = 0.0, 1 5.2, and 27.9 cm (cont) 
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Table I 

Ejector Results at x = 0.0, 1 5.2, and 27.9 cm (cont) 
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Table I 

Ejector Results at x = 0.0, 15.2, and 27.9 cm (cont) 





>r 

04 P 

fO 



OJ 

fO 

OJ 

CO 

<0 

cr Oi 


X 



X- 

p 


OJ 


X 

.4 

O P 

fn 

p 

p 

9 

o 



m4 


>0 


oi P 

a 



m 

X 

9 

X 

9 X 


m 



X 

o 

9 

P 

P 

Oj 

p 

P 

X 

N- 

p 

O P 

p 




« 

<« 

>0 


CD 


s 


P 

P> 

p 

O P 


9 



K 

p 

9 

O 

h> 

p 

p 

X 

<M 

P 

X 

OJ 


O' 



fSJ 

rsi 


OJ OI 

Pi 

OJ 

Oi 

Oj 

fO 

«i4 

fO 

9 

X p 


9 



o- 

9 

X 

X 

P 

p 

o 

P 

O 

m r- 

p p 

OJ 




fCl 

P> 

p> 


fO 

fO 

ro 

#o 

fO 


OJ 


^ t4 

o 

oj 


o 




-* 




o 

o 

t- 

o 

o o 

o 

o 


O 

O 

O 

o 

O 

o 

O O 

O 

O 

o 

o 

•o 

o o 


o 


0. 

0. 

j. 

o 

o 

O 

o 

o 

o 

1 

1 

1 

1 

1 

1 

o 



O 

-c 

IP 

OJ 

•M 


o- 

c- 

ec 

o 

o 


K» OI 


X 



o 

ft| 

sT 

p 

O' 


X 

p 

p 

OJ 


9 

f'- 




'O 

r' 


•>4 

0- 

o> 


o 

cr 

f- 

c 

>f 


a P 


X 



p 

O 

O' 


O' 

O' 

-ti 

o 

o 

ik4 


P 

9 

o 




IP 

•f 

0* 


o- 

OI 



fO 

o 


p OI 


O' 



p 

st 


OJ 

O' 

O' 

O' 

OJ 

p 

o 

p 

OJ 

OI 

OJ 



to 

rr 






p 

p 

>0 

P 

0- 

ftl 

— O 


>r 



o 

O' 

N> 

p 

OI 

0- 


P 

wm* 

X 

p 

p 

ft 

OJ 




fp 

#n 



1^ 

fO, 

fO 

fO 

fO 

OJ 

OJ 

OJ OJ O 

OJ 


C 

<tl 

OJ 

OJ 

OI 

OJ 


t- 

o 

o 

-* 

o 

O C 

o 

o 


o 

o o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

c 

0. 

0. 

o 

o 


o 

o 

o 

C 

o 

O 

r> 

o 

o 

1 

1 

1 

t 

1- 

1 

o 


'T 

4“ 

O 





c- 

P 

OI 

.-o 

o 

ft. 

— * OJ 


X 




fO 

X 

p 

P 

p 

p 

X 


O' 

9 

r- 

OJ 

o- 





Ip 


OJ 

p 

OJ 





o 

OJ 

a> -4 


X 



p 


OJ 

p 

OJ 

p 

N 

ft 

o 

O' 

9 


p 

0- 



O 


m 

P 

0“ 

p 

® 

OJ 

O' 

o 

cr 


OI O- ft 


O' 



p 

P 


p 

m 

OI 

o 

O' 

0- 


P 

OJ 

OJ 





IT 

m 

p 

p 

p 

0- 

a 

o 

Oj 



O X r- 


o- 



X 

P 

OJ 

X 

P 

O' 

OJ 

p 


Og 

P 

p 

p 

OJ 




m 


P> 


fO 

fO 

rr, 



p 

ft 

fO 

OI OJ o 

— 


o 

ft 

ft 

ft 

OJ 

Oj 


-• 

o 

o 

-* 

O 

O 

O 

o 

o 


o 

C' 

o 

O 

o 

o 

O 

O 

o 

fO 

o 

<o 

o 

o e 

o 

C7 


o 

o 

c 

n 

r> 


n 

o 

n 

1 

1 

1 

1 

1 

1 

o 


*4 


IP 



o 

fO 

fO 



OJ 



.-ti 


P 




p 

K 

O' 

.t) 

X 

OJ 

<t 

0* 

L*\ 

P 

cr 

9 

p 




•# 



o- 

0“ 

03 

O 

p> 

C 

p 

o 

p O- P 


O' 




o 

o 

•»■ 


K 

p 

o- 

X 

0- 

o 

h 

CC 

p 



er 



fp 

<o 



O 

OJ 

« 

Oj 


p 

<r fo 


(T 



m4 

st 

p 

K 

X 

X 

o 

X 


f> 

p 

ft 

ft 

IM 



ir 



>£ 

o o- 

« a 

Oj 

P 

o- 

O' 

p 

Oi 


P 



p 


o- 

ft 

9 

P 


K. 

P-4 


p 

p 

ft 

OJ 



m 



P' 



ft' 


•p 


p 

#r- 

fO. 

IT'. 



c- 


ft 

OJ 

OI 

— f 

-4 

■-4 

o 

o 

-4 

o 

o 

O 

o 

o 


• 

• 



« 

• 

« 

• 






• • 






• 

• 

• 

• 



• 


• 







0 

0 

o 

O 

0 

0 

o 

O 

o 

o 

0 

0 

o 

0 

0 

o 

O 


c 

o 

D 

o 

o 

O 

o 

c 

c 

o 

1 

1 

i 

1 

1 

o 



>0 

o 

P 


c 

p 

Q 

c 

p 

O 

X 

m 

p — 


P 



o 


OJ 

o 

o 

o 


ft 

Og 

P 

p 

p 

o 

9 



O 

cr 

nO 

K3 


r- 

c 

OJ 

p 

fO 

<> 

r«“ 

OI 

ft* fO 


P 




rr 

ft 

OJ 

0* 

O' 

p. 

o 

er 


*4. 

p 

o 

p 





1^ 

P 




9^ 

•«4 

o 

ro 

Oi 

OI p 


P 




o- 

X 

o 


p 

ft 

ft 


ft 

p 

p 

p 

p 



•C 


*0 

P 

p 

o 


9 

OJ 

p 

o- 

O 


*t OJ 


P- 




O' 

p 


X 

p 


X 

OJ 

•-4 

p 

p 

ft 

Oj 






m 

O*; 

m 


m 



■X 

P- 

ftl 

fT\ ft. 

c 



c 

ft 

OJ 

OJ 

Oi 

-• 



o 

o 


c 

C 

ro 

c 

c^ 


o 

c 

c 

r> 

o 

o 

o c 

c 

r 

o 

o 

O 

o c 

c 

o 


o 

o 

o 


O 

o 

o 

o 

c- 

c 

1 

1 

1 

1 

1 

o 




0“ 

ro 


p 

c 

o 

o* 

c 

••4 

o 

o. 

P <r 


p 



p 

ft 

O' 

c 

X. 


p 



N- 

X 


o: 

0* 



c 

c 

P 


p 

4“ 

fO 

Oi 


K 

<o 

'O o X 


o» 




o 

N. 

O' 


p 

o. 

OJ 

0- 

fO 

p 

OJ 

OJ 

p 





IP 

P 

{> 

vf 

f<\ 

<o 


OJ 

P 

o» 


O- -4 


fO 



o. 

o» 

o 

o 

O' 

X 


p*> 

OJ 

P 

X 

o 

OJ 

OJ 







•c 

o 


C“ 


p 

P 

c 


P ^ 


p 



(t 

0 


c 

K 

X 


X 

fO 

O' 

p 

p 

X 

OJ 




0Ti 


<n 

m 


ro 

fO 

>#■ 


P 

p- 


fO ft 

o 

t- 


o 

ft 

OJ 

OI 

OJ 

— 

f 

— 

o 

o e 

o 


c 

o 

c 


0. 

0. 

o 

o 

0 

0. 

o 

O 

o 

o o o 6 

o c5 

c 

o 


o 

o 

o 

0. 

0. 

c 

o 

o 

o 

o 

1 

1 

I* 

1 

1 

o 


o ir> 

o 


Oi 

p 

0y 

o* 

O) 

X 

P 

fO 

<c 

Cl ^• 


U' 



o 

o 

P 

o o 

OJ 

p 

p 

o 

P 

p 

C" 

p 

o 




<D 

CT; 


o 


« 

X 

fO 

p 

c 

OJ 

fO 

r fo 


O' 



o 



o- 



o 

cr 

p 

o 

pp^ 

ft 

s. 

O' 



>#■ 

st 

•p 


cr 


O 


cr 

c 

<T 

X 

wi4 

O CD 


9 



X. 

OJ 

o- 

o- 

o 

ff 

p 

f*.. 

P-4 

O' 

o- 

p 

X 





>0 

>0 


o o- 

® 

0^ 

•>4 



o 

X 

O- P 


OJ 



Oj 

X 

X 

O' 

p 

OI 

o 

X 

P 

P/ 

p 

X. 

f\J 

OJ 





m 

ro 

fO 

fn 

fo 

O' 

•X 


>t 

<r 

fO 

n #•*. 

o 



o 

ft. 

Oj 

OJ 


_4 



e 

o 

o 

r> 

o 

c 

o 

o 


• 

• 

• 

• 

• 

* 

« 

• 

• 

• 

• 

• 

• 

• • 

• 

• 


• 

« 

» 

• 

• 

* 

• 

* 

• 

• 

• 

• 

• 

• 

• 

• 


o 

o 

O 

O 

c 

o c 

o 

o 

o 

o o 

0 

0 

0 

c 

O 


c 

rj 

O 

o 

o 

o 

o 

o 

cr 

c 

t 

» 

t 

1 

1 

o 




-o 

Cr 


p 

C 


X 

c 

9 

't 


OJ .t 


P 




P 

-c 

X 

c 

0, 

p 

X 

3 

p 

OJ 

p 

p 

c 



o 



>C 

't 

p 


p 

K. 

p 

OJ 

X 


O p. 


0“ 



P 

o 

r 

OJ 



ft 

X 

ft. 

h» 

p 

O' 

p 

X 



tf 

«• 

IT 

P 

< 

ff 

o. 

c 

c 


ft 


>«• 

a C 


c 



r 

p 

c 

f 

X 

o, 

0 

r 

o 

c 

X 


9 

X. 




'C 

<0 

<5 

p 

p 


r 



C 


>c 

p p 


0- 





o 

r~- 

c 

K. 

p 

p 

O: 

ft 

O' 

o 


f-t 


r 

r' 

fr 

P 


r*. 

fO. 


ft- 

r" 


P 

ft 

fO 

ft fO 

p 

o 


c 

r-j 

OJ 

— 

t- 

o 

c 

c 

o 

C 

o 

r* 

r" 

c* 

C 

o 


c 

c.- 

c 

c 

c 

o 

o 

o 

c 

c 

c 

C3 


0. 

0. 

c 

c 

u 

(D 

C 

o 

o 

C 

c. 

c 


C7 

c 

1 

1 

1 

1 

t 

o 


cr 

tP 

c 


o 

OJ 


p 


o 

ft 

r 

_ 

o* ^ 


c 



f»» 

r> 

X 

p 

o 

p 

O 

Oj 

p 

P 


p 

c 



or. 


r ‘ 


>t 

p 

c- 

•4 

o 

p 

►- 


0* 

p p 


rr 




(t 

ir 

P 

[t 

r- 


O 

O' 

X 

P 

(0 

X 

O' 




m-t 


c 

9 

or 

o- 

T 


ir 

p 

V* 

p 

— < X 


c 




ft 

'T 

o, 


r , 

p 

P 

X 

r 

p 

O 

o 

r- 






>c 

P 

p 

P 

p 

p 

%c 

'O 

'X 

»o 

f»- 


p 




p 


X 



|>S • 

O' 

_ 









m 

p* 

r“. 

o , 

O') 

r* 


.t3 

ft 

ft 

ft 

ft 

r', ff* 

c 

c 


o 

o. 

— 

— 

C 

o 


o 

O 

C 

c 

c 

c 

c 

C 

c 


o 

O 

-■ 

r? 

O 

<~ 

o 

c 

o 

r- 

r- 

c 

<“ 

r- <r 

'-■ 

o 


rr 


C' 

O' 

o 

o 

r" 

r 

r*. 

c 

t 

1 

■ 

1 

1 

c 


h* 


P' 

>c 


p 


p 

o 

:* 

p 

ft 

0- 

O P 





ft 

p 

o.- 

O- 

o 

O 

p 

X 


a 

o. 

p 

p 

o 


u. 

o 

or 

•P 

cr 

og 

O" 

r*’ 

O- 


— • 

p 

p 

p 

o ^ 



U- 


p 

r 

X 

K. 

o 

P 

o. 

ft 

ft 

ft 

..4 

o 




u 

o 


0“ 

p- 

P 


o; 

o 

•M 

p 

C 

c 


tg X 



u. 


o. 

ft 

p 

P 

ft 

O 

(D 

® 

X 

o 

X 

o- 

X 



o 


If 

P 

P 

p. 

p. 

p 

p 

p 

p 

p 

p 


X Tn, 






c 

st 


o- 

P 


ft, 

o 

_4 

o 


o 

n 




fn p" 


r* 

m 

ft 

r 

ft 

ff' 

ft 

ft 

X 

X r~ 

c 

X 

c. 

OJ 


•>4 


c. 

c 

o 

o 

rr 

o 

c 

o 

o 

c 

o 


• 

• 

• 

• 

• 

• 

• 

• 

• 

* 

• 

• 

• 

• • 

• 

• 

•V. 

• 

• 

• 

• 

• 

• 

• 

• 


• 


• 

• 

• 

• 

• 



<T 

C 

O 

O 

rr 

C 


r* 

C: 



‘C 

<r c 

rr 

o 

o; 

p 

r: 

r> 

r* • 

fC 

r 

fC' 


cr 

o 

» 

\ 

» 

) 

» 

o 


















_ 



















-4 

*4 

O' 

X 

p 

p 

p 

•c 

P 

p 

P 


X 

X 

C 

^ •— • 



OJ 

ft 

p 

<► 

p 

p 

P 

p 

p 


X 

OT 

C 

X 

P 

P 

o. 


X 

X 

o- 

p 

P 

p 

X 

c- 

o- 

X 


X 

P 

p 

Oj 


X 

X 

o 

P 

P 

p 

X 

O 

o- 

X. 

,-4 

o 

O' 

X 

p 

p 

o- 

X 

e 


o 

p 

9 

o.- 

O; 

X 

P 

X. 

f*.' 

X 

P 

P 


X 

c 

..4 

Ov.' 

p 

9 

OJ 

Og 

X 

P 

1 * c 

p 

C 

'*• 

r 

p 


ir. 


p 

r' 

O 

•->4 

»-4 



1. t- 

P 

X 

P 

r 

P 

■~. 

p 

•— • 


r- 


p_ 

••4 



e 

• 

C 

• 

«4 

« 

^4 

• 

IV 

fv 

X 

X. 

P 

'* 

•P 

P 

P 

p 

p 

P *4 

C 

c 


• 

OJ 

• 

OJ 

X. 

X. 

P 

p 

p. 

P 

P 

P 

P 

• 

P 


I I I I < I I I I I I I I i < * u I t I I I t I t I I I I I I I I 


o 

C 


L- 


52 


O/.T'jq 0. 0.12^'»5 0. 1*325 0. 15565 0. 15995 0. 199HR 0.2*98fl 0.29988 



Table I 

Ejector Results at x = 0.0, 15.2, and 27.9 cm (cont) 
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Table I 

Ejector Results at x = 0.0, 1 5.2, and 27.9 cm (cont) 
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TABLE II 

SUMMARY OF COMPUTER TIME REQUIRED FOR TEST CASES 
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Table III 

Data Deck for Case I (cont) ■ 
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INITIAL TDTAL TEMPEHATUKE PROFILE 
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Data Deck for Case II (cont) 
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Table VI 

Data Deck for Case IV (cont) 
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Table VI 

Data Deck for Case IV (cont) 
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Table VI 

Data Deck for Case IV (cont) 
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Table VII 

Output from Case I (cont) 
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Table VII 

Output from Case I (cont) 
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Table VIII 

Output from Case II (cont) 
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Output from Case IV (cont) 
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Output from Case IV (cont) 
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Output from Case IV (cont) 
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Output from Case IV (cont) 
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Output from Case IV (cont) 
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Figure 1 . Schematics of Ejector in Wing Configurations 







Comparison between COMOC Predicted and Experimental Integral Parameters 
Turbulent Boundary Layer Developing in a Duct 









Velocity Profile at x equal 241.3 crn for the Turbulent 


Boundary Layer Developing in a Duct 
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7. Comparison between COMOC Predicted and Experimental Integral Parameters for the Turbulent 
Boundary Layer Developing in a Duct for Gradient and No Slip Wall Boundary Conditions 











Figure 9. Comparison Between COMOC Predicted and Experimental Centerline 
Turbulent Kinetic Energy for the Rectangular Nozzle 




Figure 10 . Comparison Between COMOC Predicted and Experimental Shape of T.K.E. 
Profile at x/D Equal 70 for the Rectangular Nozzle 
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Centerline Velocity, ft/sec 



Figure 11 . Comparison Between COMOC Predicted and Experimental Decay of Centerline 

Velocity for the Axisymmetric Nozzle 
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Exit Plane (x = 0) 


Figure 12. Hypermixing Nozzle Geometry 


u (ft/sec) 



x = 5 in.; V = 0 in. x = 5 in.; y = 0.4 in. 

(Mid Element) (Between Elements) 


Figure 13. Comparison with Eastlake Hypermixing Data 
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Figure 14 . Configuration and Definition of Parameters for the Wall 

Jet Developing in a Duct 
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F igur e 1 5 . Comparison Between COMOC Predicted and Experimental Skin Friction and Velocity 

Profile The Wall Jet Developing in a Duct 










(a) Maximum TKE versus x 
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Fi jure 17 . Comparison Between COMOC Predicted and Experimental TKE Profile for 

The Wall Jet Developing in a Duct 
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(a) 



(b) 

Figure 19. Developing Flow in a Pipe 
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x/a=-20 


x/a = 30 


Figure 22. Axial Velocity Profiles for Conical Center Pipe Flow 






Figure 23. Maximum Secondary Velocity Comparison for Conical Center Pipe Flow. Discrepancies Between 
The Two Codes In the Initial Region are Due to Differences in Descretization and Axial Integration Stepsize 
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Figure 25. Growth of Centerline Velocity for Fluted Pipe Flow 
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Figure 29. Axisymmetric Fluted Mixer Geometry 
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Figure 3 1 . Static Temperature Profiles for Fluted Mixer 





